{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\erihe\\Anaconda3\\lib\\site-packages\\ipykernel_launcher.py:159: RuntimeWarning: invalid value encountered in true_divide\n"
     ]
    }
   ],
   "source": [
    "import * from utils\n",
    "_SQRT2 = np.sqrt(2.0)\n",
    "_SQRT3 = np.sqrt(3.0)\n",
    "\n",
    "npoints = 1000\n",
    "sigs = [0, 1, 2, 3, 4, 5, 6,7,8,9,10]\n",
    "num_peaks_all = {}\n",
    "cs = np.array(['w','g','r','c','m','y','b','k',])[:, np.newaxis].repeat(5,1).T.flatten()\n",
    "\n",
    "numbins = 50\n",
    "bins = np.linspace(0,2*np.pi, numbins+1)\n",
    "numangsint = 51\n",
    "numangsint_1 = numangsint-1\n",
    "mid = int((numangsint_1)/2)\n",
    "indstemp1 = np.zeros((numangsint_1,numangsint_1), dtype=int)\n",
    "indstemp1[indstemp1==0] = np.arange((numangsint_1)**2)\n",
    "indstemp1temp = indstemp1.copy()\n",
    "for rat_name, mod_name, sess_name, day_name in (('R', '1', 'OF', 'day1'),\n",
    "                                            ('R', '2', 'OF', 'day1'),\n",
    "                                            ('R', '3', 'OF', 'day1'),\n",
    "                                            ('R', '1', 'WW', 'day1'),\n",
    "                                            ('R', '2', 'WW', 'day1'),\n",
    "                                            ('R', '3', 'WW', 'day1'),\n",
    "                                            ('R', '1', 'OF', 'day2'),\n",
    "                                            ('R', '2', 'OF', 'day2'),\n",
    "                                            ('R', '3', 'OF', 'day2'),\n",
    "                                            ('R', '1', 'REM', 'day2'),\n",
    "                                            ('R', '2', 'REM', 'day2'),\n",
    "                                            ('R', '3', 'REM', 'day2'),\n",
    "                                            ('R', '2', 'SWS', 'day2'),\n",
    "                                            ('R', '3', 'SWS', 'day2'),\n",
    "                                            ('Q', '1', 'OF', ''),\n",
    "                                            ('Q', '2', 'OF', ''),\n",
    "                                            ('Q', '1', 'WW', ''),\n",
    "                                            ('Q', '2', 'WW', ''),\n",
    "                                            ('Q', '1', 'REM', ''),\n",
    "                                            ('Q', '2', 'REM', ''),\n",
    "                                            ('Q', '1', 'SWS', ''),\n",
    "                                            ('Q', '2', 'SWS', ''),\n",
    "                                            ('S', '1', 'OF', ''),\n",
    "                                            ('S', '1', 'WW', ''),\n",
    "                                            ):\n",
    "        \n",
    "    file_name = rat_name + '_' + mod_name + '_' + sess_name + '_OF'\n",
    "    if len(day_name)>0:\n",
    "        file_name += '_' + day_name\n",
    "    \n",
    "    f = np.load('Toroidal_topology_grid_cell_data/Results/' + file_name + '_toroidal_rate_maps.npz', allow_pickle = True)\n",
    "    mtots = f['mtots']\n",
    "    f.close()\n",
    "    num_neurons = np.shape(mtots)[0]\n",
    "\n",
    "    num_peaks = np.zeros((num_neurons, len(sigs)))\n",
    "    for nn in np.arange(num_neurons):\n",
    "        mtemp1_3 = mtots[nn,:,:].copy()\n",
    "        for j in range(numangsint_1):\n",
    "            mtemp1_3[j,:] = np.roll(mtemp1_3[j,:],int(j/2))\n",
    "        mtemp1_4 = np.concatenate((mtemp1_3, mtemp1_3, mtemp1_3),1)\n",
    "        mtemp1_5 = np.zeros_like(mtemp1_4)\n",
    "        mtemp1_5[:, :mid] = mtemp1_4[:, (numangsint_1)*3-mid:]  \n",
    "        mtemp1_5[:, mid:] = mtemp1_4[:,:(numangsint_1)*3-mid]  \n",
    "        mtot = np.concatenate((mtemp1_5,mtemp1_4,mtemp1_5))\n",
    "        for s, sig in enumerate(sigs):\n",
    "            data = gaussian_filter(mtot, sigma = sig, mode = 'constant')\n",
    "            datatmp = data.copy().flatten()\n",
    "            datatmpBU = datatmp.copy()\n",
    "            datainds = np.indices(np.shape(data)).flatten().reshape((2,data.shape[0]**2)).T\n",
    "            datacount = np.zeros(len(datatmp))\n",
    "            XY = np.zeros((npoints,2))\n",
    "            for k in range(npoints):\n",
    "                ktemp = np.argmax(datatmp)\n",
    "                XY[k,:] = datainds[ktemp,:]\n",
    "                datacount[ktemp] += 1\n",
    "                datatmp[ktemp] = datatmpBU[ktemp]/datacount[ktemp]\n",
    "\n",
    "            X = squareform(pdist(XY, 'euclidean'))\n",
    "            thresh = 5\n",
    "            X[X>thresh] = -1\n",
    "            knn_indices = []\n",
    "            knn_dists  =[]\n",
    "            F = np.zeros(npoints)\n",
    "            for i in range(npoints):\n",
    "                indsvals = np.where(X[i,:]>=0)[0]\n",
    "                knn_indices.append(indsvals)\n",
    "                knn_dists.append(X[i,indsvals])\n",
    "                F[i] = np.sum(np.exp(knn_dists[i]))\n",
    "            i = np.argmax(F)\n",
    "            inds_all = np.arange(len(XY[:,0]))\n",
    "            classes = np.zeros(len(XY[:,0]))\n",
    "            classcurr = 1\n",
    "            inds_left = inds_all>-1\n",
    "            inds_left[i] = False\n",
    "#            inds_left[knn_indices[i]] = False\n",
    "            classes[i] = classcurr\n",
    "            classes[knn_indices[i]] = classcurr\n",
    "            for j in range(npoints-1):\n",
    "                F[knn_indices[i]] -= knn_dists[i]\n",
    "                Fmax = np.argmax(F[inds_left])\n",
    "                i = inds_all[inds_left][Fmax]\n",
    "                inds_left[i] = False\n",
    "#                inds_left[knn_indices[i]] = False\n",
    "                if classes[i] == 0:\n",
    "                    classcurr += 1\n",
    "                    classes[i] = classcurr\n",
    "                    classes[knn_indices[i]] = classcurr\n",
    "                else:\n",
    "                    classes[knn_indices[i]] = classes[i]\n",
    "            ind = classes.astype(int) -1\n",
    "            binstacked = np.bincount(ind)\n",
    "            indd = np.unique(ind)\n",
    "                        \n",
    "            inds1 = ((XY[:,0]>=numangsint_1) & \n",
    "                    (XY[:,0]<2*numangsint_1) & \n",
    "                    (XY[:,1]>=numangsint_1) & \n",
    "                    (XY[:,1]<2*numangsint_1))\n",
    "\n",
    "            inds2 = ((XY[:,0]>=1/2*numangsint_1) & \n",
    "                    (XY[:,0]<3/2*numangsint_1) & \n",
    "                    (XY[:,1]>=numangsint_1) & \n",
    "                    (XY[:,1]<2*numangsint_1))\n",
    "\n",
    "            inds3 = ((XY[:,0]>=3/2*numangsint_1) & \n",
    "                    (XY[:,0]<5/2*numangsint_1) & \n",
    "                    (XY[:,1]>=numangsint_1) & \n",
    "                    (XY[:,1]<2*numangsint_1))\n",
    "\n",
    "            inds4 = ((XY[:,0]>=numangsint_1) & \n",
    "                    (XY[:,0]<2*numangsint_1) & \n",
    "                    (XY[:,1]>=1/2*numangsint_1) & \n",
    "                    (XY[:,1]<3/2*numangsint_1))\n",
    "\n",
    "            inds5 = ((XY[:,0]>=numangsint_1) & \n",
    "                    (XY[:,0]<2*numangsint_1) & \n",
    "                    (XY[:,1]>=3/2*numangsint_1) & \n",
    "                    (XY[:,1]<5/2*numangsint_1))\n",
    "            nump = []\n",
    "            for inds in [inds1,inds2,inds3,inds4,inds5]:\n",
    "                bintemp = np.zeros(max(ind)+1)\n",
    "                ind1 = ind[inds]\n",
    "                inddd = np.unique(ind1)\n",
    "                bintemp[inddd] = np.bincount(ind1)[inddd]\n",
    "                indsfinal = np.where(np.divide(bintemp, binstacked)>0.5)[0]\n",
    "                nump.append(len(np.unique(indsfinal)))\n",
    "            \n",
    "            num_peaks[nn, s] = np.max(nump)\n",
    "    num_peaks_all[rat_name + '_' + mod_name + '_' + sess_name] = num_peaks\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x168af056608>]"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAOwAAAD4CAYAAADvlAqZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzdeXhc1Xn48e+5s89o12jfbXmTF2wjNkPAYPYQyNIEHCAhSUuaJs3arG1omrTNr20S2qYJCU1SEkgDBKghxMR4A2xjvNuSvFuyJGuf0Tr7cu/5/TGysIzBAo89M+J8nuc8o5k5c+87sl6fc88991whpURRlMygpToARVGmTiWsomQQlbCKkkFUwipKBlEJqygZxJyqHbvdbllbW5uq3StK2tq1a5dXSll0pvdSlrC1tbXs3LkzVbtXlLQlhOh4s/dUl1hRMohKWEXJICphFSWDqIRVlAyiElZRMshZE1YI8SshxIAQouVN3hdCiP8UQhwTQjQJIZYmP0xFUWBqLewjwM1v8f4twKzxcj/w0LmHpSjKmZz1PKyU8hUhRO1bVLkD+I1MXKf3mhAiTwhRJqXsTVKM04eUxA+9TGT1QxALT/VDYMRBjxONxBj1wWAwj6A+B4t1JkK881PpybqwUgoBQpy2TfG2ozGkjkEcQ8aRGKds4e1FKicV+Q5iOQMBwqRhRUzamphCbPK0J7X3XkblZY3vKIxkTJyoAE6c8rxr/LU3JKwQ4n4SrTDV1dVJ2HWGGDgELU9h7HgMc6gXM3D6ZcjSgHhYIx40EQuaiAQtjETLGTaqGKEOv7UKW14dea4CSvNNlGjjCaKuZ844g0ePpTRhz/Tf1xn/iqSUDwMPAzQ2NmbUX1pwwMPWR1cRj8XOXjkexRrooXB0P+XHmrDpAWKldoZCWQxE5jBmcxLVLcggWEI6FsOB2ZRLyFlCwFVGwFlGMLcYq8lCiUVQataYaxGYhCAm4wyaTtCXu4+mqtcYzfKN71QgsWBgIgaEdYOgHiciIWZATIIwrNilCztObIYDS9yGOWaGqMRv6ATNJkIOG0G7A7/Vjt/sICCcyFP+iU3EyZZ+so0Q2fEoWREdZwgsfjMmvx2h2xCGgSUawRINY42GsUZDWCNhLNEImmFMbEsChgmkGbCb0FxWzLl2LMU5eKwx2mLD9As/o+YgYXMAQ0Qm/ZolGoapAMNUgCZysRrZOHQnrpiNvLCVwpBGdkzHYophNkcxm6OYzFHMpgiaKYrFFMFsjmI1RTCbI5jNMUzWGGZLDM1sTNpX5AToWzWqm0OUDcbQBTTVCbY0CHbW5WO35jLDpnNxzhBzsgbR3uJgc+bcP5z9b+hNJCNhu4CqU55XAj1J2G5akLEYw7/7HSce/E/KQ4G3/XlH5yAA/dUlADitgrGSRYSzKvGVVRC2F4H2+j9DjmWMakeMEkuUXN0JQNA6wuH8Q+y3HeJAtI3gSJwsj0Z2hxXNKJy0Pw2wjZfcM38jIDBeTonzbX+zM3OGzWSFTJM6jgF7nEFXjLHiOGOuGKOuGKNZcQKOOPIteqtWzYoz7iQ36qAsmI9Dd5Add5EjcrG5ypDZlfhc2QxaHAyarAyaLPSZzIQ0U5K+zSmqIbd0jPKbPTR2H+HSvfuobzrKkueDRM1D7Jo5ypb5Br+aKYgNFZMnapmdU0uju5LFednEYhGCgSCBQISSwnfeu0xGwj4HfE4I8ThwGTA6XY5f/Zu30P/97xNtbaWleDadNzfyCfdrMHICIqMA6AJMljgOIXFGQwAEHA568rPoyrNQ9C+J/2oPft1NoPtWfF0XIw0zmssHucPEs1qwyziFkSxmBisoiicSsF12soMddAWPEh0YIvuIGU0K5uICIGizELXaJrVYCQJNmDBpZjRhQqChoSUOMYUEYSAxMISBIXQMYQAGmjTQdANNB2loSPnGJkLTNIQQaJqGZtLQhIZJ0xCalnhN07AUu7Dn52M7pZis1rf1e5fhMFnWLArcZdjtdhwOx8Sj2Xz2P1lfXKczHKU9FGEkpr+tfRuGQSwWm1SisSjecJRDuoU2VxWPNszk0YZbYKWk8VAzN217hcua93LF4QAhq5nd9SZennecnfXNbPcLpGHBRR0zs+dzZeUlRAzxjv+DFGc7BhJC/A5YDriBfuDvAQuAlPJnQggB/BeJkeQg8Akp5Vln9Tc2Nsp0nfwf7eig///9C/6NG7FUV/OHKz7AzyIlPFX5ANnlLkJ2C1bfMLmeQQoGg2gSgnaN/iIbw5UVGLl1yL09yP1eGp4fJWgp4skbHkQagiPuHRwt2sSSUCWX+y+iITYbB3biMkZvuJ2+QCs9wWOE9QC6Jgk6JSGHlb78CrrL5hOyFzC3r4ebghFKXaXYRS7WqANLzIotZsYmBUK83myNGGF6LX6GLH5GzH5G5BghPTTxfmFhIWVlZeTk5ExKjNMfbTYbJtN5aLky0GgsTrMvxI7BEfaM+NgfjNIb01lyZD8rdrzK1Xu24wqH8Dkc7J5RyqZZdlpqfOiOAYQw+M3yZ1hSM+tNty+E2CWlPONB7lRGiVee5X0JfPZs28kEuj/A4M8eYvDXv0GzWMj+3F3svyiH//xjBbfWrkN3j2Dr66dyMIZZN4jbnfjnXUV83i1YaldQ7ayh4IWN9P71dwhh47Xln6E8+j8QERwr2IO9NsD17XP4avcVaMJEKO6nJ3SYtthhDpmPMuaKYnO6yDPyKDfqaa+4mOdnVTHgtFDn1/n48Sg3HY1jkRWJgBONPFFDEtTAa9cZcQTx2QOMCh9DoUGCoeDE93PnuZlVPouysjLKysooLS3Fbren4Ded2XItZq4qyOaqguyJ18biOi0Xz6Hptpv59dAocstm5m95mSuad3PN/hgDBYXsWXQ1u2sLKHdXvON9n7WFPV/SqYWVhsHoqmcZ+NEP0L1DGNeWM3TLMNGsUR7c9Zd0jFTxcs0PKOg5irTnIubdDgs/DLVXwfjxUnv3fk585wFyXutg77wbGCm8Gk1aeP+Ln8dqN9Hy4f9iTkQjZkQ4ENnDC+5t9DtHqBI1NERn0RCcQWW0hGGrxhM1Zp6qsjFm0Vg0qvP+YxHmdUeJGxAXYC+04yp1EMwew2/2MRwcpK+/j0AgcVwqhMDtdlNeXj4pOW02Wyp/ze86/rhOS7+X3hfXYl/7IlX7dmPSdUpXryZ/Rt2bfu6tWth3dcJKqePd8gyD//pj5BEP0VqD0Y/oiNmFFBZcw/Emnc/uvZbni3/OgrFX4Kovw/JvgNlG3Iizd2AvL3e9TOfLL/D+p0bw5a+gs3o5YKV4oZXy/HwKXtlHgQ4jbjcHxrZyzLeXUFUBc/0VFOeV4Cp143AX0W0r4HHNxHqHTkzAnO4Yyw6FmStNlNblUlKXQ+mMXHJLbTQ172Pr1q2MjIwghKCoqOgNyWl9m8eNyvkXHx5meMurFN323respxL2FNHoEENDmxg8uoboLzdhfy2OniuJr6wk53134C5aTnb2fNj+S+5YFeTD1q3cK16AKz6Hb/nX2dL7Ki+deInN3ZvxB0dYucnJ/L7lnKi6FkOzMntxMfkyTl6nH7uEMc3Po+7nMY4exT1swVlawpW3/jm6Uc5Au49dgz42VJo4WGlFk3DZsOQucxaNtfmU1OXgyku0isFgkO3bt7Nt2zZCoRCVlZVceeWV1NfXY7FYLvjvUTl/zukYNtNJKfH5mvEOvsTg4MuMDe4ja70ge40Zm65hv/tayv/6W9jyKl//0OEX+OMfn+IarZJ7xQvsqV7Cv8XaOPjENcRlnDxbHrfrlzJ3bT49rqV01NiZPSeX0ohObvswGoIxLcDPyp4l3LaaizbkMuaqxpl/EUbkOjY/E+Z4SSfbF7k4OseJE8En8nL57JwyKpyTu63Dw8Ns3bqVPXv2EIvFmD17NlddddW7a+KJMmHaJ2xX96McOfIPIAX5R+sof7wQ+sbIWrGCkq9/Dev4H37MiHF46DB7jz7Pnj2PUMhlfMvyFM+6XPyTxc98s51PLPgEVxQsI/DbLg4dt3Ai18GsfKg1W3H2J0ZeB7UIj5W9wJ9y1pIVtPDbJ3xYdD//9947KGj4IEfq7azNMjiixyixmvl2VTH3lheSY548AtvX18eWLVtoaWlBCMGiRYtYtmwZxcXFF/x3qKSPaZ+wvb1PkztSj/vZIkLbdmGtn0nprx4ksnQuWwb2sXfX0+zz7GO/dz9hPTG/9+6xMr6hraKl8GLq73yQLe4GhK6x85mD7N7QRVzkMNcRZYbTjEUKjKhBr4zzu5z1bKj4I8KQLD2cx8UD9Vj0DszOLDq/+in+wztGQI8wx27nR1WlfKgkH9spU2KklLS3t7N582ZaW1uxWq1cfvnlXH755eTmnnkahPLuMq0TNhzuIbplP4UPWwm5Bun5i1vY2Ghjb/8/0/FEYp0rs2amoaCBD8+8ncX7VrFkcICiWCe7bJey9DPPI8yJLuraBzdx9GiMBruNOrsJk7AQl5LWiM7q3N1sqn4Gn2mUGd0uruqt4bIb72LnM4+jW6x05RTytGeU24vzuae8kMYc56RzpYZhcPDgQbZs2UJPTw8ul4sVK1bQ2NiIw5GsOUjKdDCtE9bjWYdjp8aYE770yTA+51oK+gpYXLSYD836EIuLFzOvYB52BDz2Ieg/jiElm/X5OD/ym4lkPfDMTvpaY9ySY8aiCXyGwTYh8S/1s9n2GPtHWygas3P1gQrql99Hy00LWfPfP8ChG/jcpdQ6HexdNp9cy+RfdywWY+/evbz66qsMDw9TUFDAbbfdxkUXXaQGkpQzmrYJG46HefXYQ8xp1+itd/PNG7/K4uLFVGZVTmrdkBL+79PQvgkpTOyRs3mi/l/48cwyADo2NPHq2lGWOs3oQvI/bo1LbiqhdfS3rGpdhdNn4dIjpTjmfpRX72nk1+EYf/bEr8gfHWTpl/6OBV/7WmJW7SnJGgqF2LFjB9u2bSMQCFBeXs4NN9zA3Llz0d5q1rjyrjctE3YgOMDfbPgs9+oD2LxWLv/YPbhnvu/MlTf+EzQ9AZqJHvssPjH8Nzxzy2IADr7Qwr7nvZiEoMgsaKq2U3t1C//Y9BUisTA1/ZVES+5k7S0XEUWwyGblgZZNBDsOs+JTf8XiSy+FRx+dtLujR4/y+9//nmg0Sn19PVdeeSW1tbWT/xNRlDcx7RK2xdvCFzZ8gTqTF0dPIgkcixadufLuR+GVfwPNTDR/Fnf0fZlbG+cww53F5kf34dk+xKAOC51gCPjfgh/StPsA2eEKYoUfY2dNA1kmjbtKEsem5j2v8cKG1Vx0w60svvHWxD6qXr+Q6ciRIzzxxBMUFRXx/ve/n9LS0vP961CmmWmVsC8cf4Fvb/k2hfZCPlm3FOv23SBi2OfPf2PlY+vhuc+DMEF+Ld/N/yd8/VE+s6yO5x/cSey4D58BdiEptVpoKeqkOXoMX+GX8biWcJHDwt/WlHJHUR4us4meI4d48uEfUzV/Edfed//r+3niCQAOL17Mk08+SXFxMR/72MfUYJLyjkyLhDWkwU/2/oSHmx5mafFSfnjNv9Ky40byu/OxzsjClJ09+QN9zfD4PYllTXIqOHTjb3nsf1q57+IqXvmvJizDIXLMGn1xCOSCE8F/u54k7LyK9xUu4TML59CQ9XrC+Qa9PPuDfySroJD3fekbmE69BOyhhwgGgzxx222UlpZy7733qmRV3rGMH+EIxoJ8+aUv83DTw3xw1gf5xY2/QISOoMf9iNYgjoULJ39gtBt+8wHQw+Byw31/4J+3jJFtNVGweRDT4ChzrAbHo9CUB9c6bezLHqbVcZy/NWbz4ysWT0rWWCTMqn/7HvFohA987QEc2TmTdhcIBhnweCgrK1PJqpyzjG5he/w9fH7D5zk6cpSvX/J17p53N0IIPN51mEfsyGEf9oULXv9AeAwefT8EveDIg/v+yBavi1eOeFgeMlMaHqTRamJd3M6IXRBYkEVZS5jfVDxLldfFR+6/d9L+pZSseeg/GGhv4wNfe4DCysnTBQ8ePIjD48FqtXLvvfeqS9mUc5axLezegb2s/ONKevw9/HTFT7mn4R6EEEgp8XrXUTA4FzhlwEmPwe/uBO8RsLrg488TcdbyzUf3kG0Ibo90M88c4rfuQmJx+OMlTr4+KAg7Y7ySvYPbi296Q+u57ZknOLx1E+9Z+XFmLL1k0nsHDhzg97//PTarldKSEpWsSlJkZMKuOraKT675JFmWLB5772NcWXHlxHs+XzORSB+u7gKExYJtzpzxc61/CR2vgtkOH3+OETGD73z/VTqjMe6K9dNV7uTzV84hvzfO3lITv1xYS35viOdc68kOmrn7ts9PiuHo9lfZ8uRjNLznWi65/UOT3tu/fz+///3vqaiooKSkRJ1bVZImo7rEuqHz4K4H+fWBX3NZ2WX88JofkmubPMfW410HaIjWALZ589CsVlj/D9DyVGKxs3tX0T5Uw59+uZ01thClphjrL6niSFEhf/WnUcZMkrs+OJuSA8P4LZIn8l/kVu0KsvNfX+xsoL2N1f/1Q8rq53DD/X896RxqS0sLTz/9NFVVVdx9991od9xxoX49yrtAxiSsL+rj6698nU3dm1g5dyVfveSrWLQ3Tt/zeNaSl9NI5MBh8j7wAdi/Cjb9CISG/OhT7DpQxrbnmthfZmIoaBC7qIS8XBtf3hbEFjDQ31PM8poCeh9vZWvuHmJGlE+/9+sT2w+OjrDq376H3ZXF7X/zt5hPuVC8ubmZZ555ZiJZbTYbqFUelCTKiL5a51gn96y+h609W/n25d/mW5d964zJGgp1EggcoTCwGBkMJgacNj8IQOwDj7Bmg5utz7XRdnUBq+NBjDwrHwkIfnbQwNIVxlNs4XMrF+B/rRdpSH6Z938sE/MpLUkMJunxGM/96J8JjY3x/q9+m6z8gol9NzU18cwzz1BdXf16sgI88kiiKEoSpH0Lu713O19++csA/PyGn3Np2aVvWtfjWQeAs6eAEOMDTk92oJuyePrZUvZFhtn04SKOtw1jiRr8vWbmhhEzj7WNYjEJPvuFRogbBLb10prdQa9tkAev+RGQGBFe94uH6D50gPd+4WuUzKif2O++fftYtWoVNTU1fPSjH528PMvJZL3vvmT+WpR3qbRO2CcPP8n3t32fmpwafnzdj6nKqXrL+h7vOrJcc9APd6NlZWGtqUEP+zkau5Rfzxbsqc6hKBTEeWSYK4SFW+wuHhoJUBoXLFo5i6JCB/5tvRjBOL8q/z/m6VXMr0vcjG//y+tp2fgil3/wTuYuu3pin3v37mXVqlXU1dWxcuVKtZaScl6lbZf4hzt/yPde+x5XlF/BY7c+dtZkjUaHGBnZgbvoesJNzdgXLkAMtWEixqqZ82mpsfIXkVGufuwFdGHiL23ZPFNkpnggjmN2Du+5pgppSPybuxmyj7An5wh/fulnAAj5fbzy2K8on9PAsg/fPbFPlazKhZa2CbvQvZD75t/Hj6/7MVnWrLPWHxzcCBi4s64mfOQIjoWL4NiLABx1zuT5+DDXPfAd1lQ18l7NSvjKcry7hpA2Ex/99EUARI4OE/eEeMLxPMVGLjcsuA2ALY8/StjvZ8Un/xIxfopmz549rFq1ihkzZryxG6wo50naJuyNtTfylcavYJrifVI83nXYbKVYukwQj+NYtJBY+zYAbEYlti99kf99z2cxCcE9y2fymzVtFBkat9zXgN2VGMDybekhYo6yumTbxKyp/rZj7Fv3Aotvfi/FtTMA2L17N88++ywzZ85k5cqV6mJz5YJJ62PYqdL1MIODmygr+xDhbc0A2BcupP93J7Bb8qjdvJsTy/6ctRYH99W4+X87O7kuqFG9tIhZSxKLmsX6A0SODPOiYyNmqXFX48eRhsH6Xz6EMyd3oiu8c+dOnn/+eerr67nzzjvPnqyrV5/X7668u6RtC/t2DA1vwTBCFLmvJ9Tcgrm4GEtJCTLk4YBrJtcNRfnv/FnkmDS6rLCoR8fqsnDD3XMntuHf0oMhDH5bvo73Vd2K0+Kk5eV19B47zNV3fwK7K2siWWfNmjW1ZAVwOhNFUZJgWiSs17MOkymL/PzLCDc1YV+0EPQ4xdE+urU6TtRcy3Z0LqkvJNA8QrGuccO98ya6wnogRmB3P7tMexizBPjU5Z8h5Pex6bePUDG3gYarr+PQoUNvP1kBfvrTRFGUJMj4hJVSx+Ndj7twOXIsSLSjA8fCRXjbt2GTMUS4nJ+JKIVOK00HvFwZsTCrsZgZi4smthHY1gtxya/LXuCKgkuoyKqYGGi67hN/ia7rrFmzhuLiYu68884p3fJwwpNPJoqiJEHGJ+zo6B5isUHcRdcTatkPgGPRQg4ffAUA32g5hzGIxnTuiNpxuCy8567ZE5+XcQP/1h6OizZas3r480s/84aBpu3btzM8PMyNN9749pJVUZIs4xPW412HEBbchcsJNzcBYJ8/n+HeI0SFmZZYFQ5gcVAjPwLLPzoHR9brp2CCzV4MX4wnCtZQa6/i4qKlkwaaAoEAL7/8MvX19dTX179JFIpyYUwpYYUQNwshDgshjgkhvnGG96uFEBuFEHuEEE1CiFuTH+qZeb3ryM+7DLM5m1BzC9a6OmRWNtnBE3SbqzimmSnTNS4PWai/uJiZS1+/1YWUEt+mLoYZ4uWiZu5b8in2v7J+0kDTyy+/TDQa5cYbb7xQX0lR3tRZE1YIYQJ+AtwCNAArhRANp1X7O+BJKeUS4C7ggoyyBAKtBIPHKSq6ASkloaYmHIsWsnMswJxAG6PxGo5LnWtCFuxOM1ef0hUGiLaPEe8J8EfnerLNWawovppNv32E8jmJgSav18vOnTtZunSpuqeNkhamckB2KXBMStkGIIR4HLgDOHBKHQmcXI4hF+hJZpBvxuNZC4DbvYJ4Xx+614t94SK2drdzeXSQnkAJ9XET7rjG1XfNxpE9eTaSb1MXERnm6crN3D3vY+x8+qnXZzQJwdq1azGbzVx77bXvPMiXXjqHb6gok02lS1wBnDjledf4a6f6DnCPEKILWA389Zk2JIS4XwixUwix0+PxvINwJ/N415GdvQC7vYxQU2LChGPRQrqPJ+47269XURnXiFsF9RdPbiHjgyFCBwbZbHqVqEnnBvvl7Fu7emKg6fjx4xw+fJj3vOc9ZGWdfWqkolwIU0nYMy1Jf/pdoFcCj0gpK4FbgUeFEG/YtpTyYSllo5Sysaio6PS335ZIxMPY2F6K3DcAJAacLBZG6maSNXQYgAOmWRQaGo4i+xtW1ve92o2Ukt/WbOD66hXs+9/fTww0GYbBiy++SG5uLpdffvk5xckPfpAoipIEU0nYLuDUS2UqeWOX91PAkwBSyq2AHXAnI8A34/WuByRFRYmEDTW3YJ8zh43+MPMDx/DJAg5o+bh1QVn15MXTjHAc/7ZeWmQTvfYhrhtrmDTQ1NTURG9vLytWrDj3ecLPP58oipIEU0nYHcAsIUSdEMJKYlDpudPqdAIrAIQQ80gk7Ln3ed+Cx7sWu70Kl2s2UtcJt7TgWLSQDYNjLPQdwx+tos8wMCOonZE36bOBHX2IODxRsZEFrrl0/GHDxEBTNBpl/fr1lJeXs2DBgjfZu6KkxlkTVkoZBz4HrAEOkhgN3i+E+K4Q4vbxal8B/kIIsQ/4HXCflPL0bnPSxOMBhodfpajoeoQQRI8fxwgEsC5YyNbBIWYGOwlHSonpifrFla+v/C8NychL7XTHO9idd5QVnXWTBpq2bt2Kz+fjpptuUqsdKmlnStN2pJSrSQwmnfraA6f8fAC48vTPnS9DQ5swjOjE8evJAadjdfWUdB3GQpxBKikwEgmXX/b65PvQfi8iIHnevYnaYCGj2w+y5ObbKK6dgc/nY/PmzcybN4+ampoL9XUUZcoysgnxeNZiNueRm3sxAKHmJjSXi3VZ+Sz0HwOg1TQbty7QHRpW++v/Lw2uPYpPH+EPJa9x7ZGKSZfObdiwAV3XueGGG5IXrMORKIqSBBk3MdYwYngHN+J2X4emJcIPN7dgX7CADSMBPj16hLi00Gyqwa1rZFW83rpGu3yIAZ2XnVup73Ehe0a5+q++hN2VRV9fH3v27OGKK66goKDgzXb/9r3wQvK2pbzrZVwLOzK6k3h8dKI7bEQihA8fxmiYT4s/xMKhw4zGKmmVGgWGoLzq9ePXgRcOEjUiPF3xCpcdKZoYaJJS8uKLL+JwOLj66qvfbNeKknIZl7Bezzo0zUZh4XsAiBw6BLEYB2tngpRUR9oJR0vxGgYmBDXjI8T6WASjNcQe0x5qjprQIsbEQNPRo0dpa2vjmmuuSf7d5b73vURRlCTIqISVUuLxrqUg/0pMpkRX9+SA09rSKhbIEZzCT9goQ4yPELsrErOU+l84gJCCNcWvMbczmyU3JQaadF3nxRdfpKCggMbGxuQHvX59oihKEmRUwvr9hwiHu3EXXT/xWqi5CZPbzWph43ZfGwA9lnqKdA0J5Jc6MaI60b3DtBpHyD80hjnLwbKPJAaadu/ejdfr5YYbblDXuippL6MSNnGjK4HbvWLitXBzC+F5DYzqBks6dgNwQJuD29CQLhNmqwnPhkOYpYVtWbsoHrFx3T1/gd2VRTgcZuPGjdTU1DB37tw32auipI+MSlivdy25uUuwWROzHvWxMaLHj3O0rh6TgJqBg/jjhRyUObh1QW6ZK3HN6+ZuvPE+4sfa0SrzWbg8MWC1adMmgsEgN9100xvmGitKOsqYhA2He/D59lPkfr07HG5pAWBjaTWX5LjIkT0EYyUclwZ5hqCiJofBba3Y4w6atb1Yoxrvvf9LCCEYHh7mtddeY9GiRZSXl5+/wAsLE0VRkiBjDtoS3WFwu1+f1HBywOlPxRX8tcUgx9JPf3QOI8JAQ1BVm4Pn5V3YDBuejiYii9zMnpO4V8769esRQrBixYo37iyZnn76/G5feVfJmBbW61mH0zkTl2vGxGuh5mbClZX4XVlc0rwFTRiMmWtx6omvVViRhRiGwXgfQWuc2z6WuEy3q6uLlpYWli1bRm5u7hn3pyjpKCMSNhYbY3hk26TuMEC4uZnjdbMosZopOrQDgDbzPNy6QArILYKp8vwAABuASURBVLDjjLsYCfYycGk2i6sakVKyZs0asrKyuPLKCzD9+ZvfTBRFSYKMSNjBwZeQMk7RKadzYv39xAcG2FxRw7UFOWiBdmKGlX3xWty6hsixEOkaRhMm+vUebnnvJwA4cOAAJ06c4Nprr339psvn09atiaIoSZARCevxrsVqdZOTs3jitVBTYknTPVUzuDbHgcvswR8voRVBkSHIL3Mxsj+xsk2nuYvra68nHo+zdu1aiouLWbJkSUq+i6Kci7RPWMOIMDj4Cm73Ck5ddSbc1IxhMnG8uoaLO46TZ+shrrvpkAa5hkZlbQ7h9hHCepCxnCgWzcL27dsZGRlR17oqGSvt/2qHh19D1/0Tk/1PCjU301Vdy0WF+YS3bMGuBZBaMQE9cd18ZU0ucjDOUKSXrNKiSQuCz5w5MxVfRVHOWdonrMe7DpPJSX7+sonXpGEQamlhT2Ud1xXkEO5MrLjaZ2ogb3yEOL/Iji1qYyjSR1F5TeoWBK+sTBRFSYK0Pg8rpYHXs56CgqsxmV4fIIq2tyP9fg7VzuQr+S5EvB+AJmMRRYbA0MAeiRNAYyjSS6F7CTu37OTiiy++8AuCP/bYhd2fMq2ldcL6fC1Eov1vOJ1zcsDJM2sO9d2d9LvGCMTzOWxk49bBkmcl0JpYA2442ofo8mE2m1m+fPmF/gqKklRp3SX2eF5ECBNu9+SV94NNzYRsdmY3zMW3Yw+59n5iRiGtGLh1QUF5FsFWLz59lCHzGP3HvalbEPyLX0wURUmC9E5Y7zpycxuxWCYvUzq0dy+HamZwbXE+/bsOkWfqRQo3JwydbKlRXZeL0R/BE+8l7BLYbfZzXxD8ndq7N1EUJQnSNmGDwXYCgaMTC4WfZESjcOQoR2pmcHVeFmFPB0JIwrIK3UhccVNa4sQcNjMc7iXmsFBRUXHuC4IrShpI24Q9Odn/9NM5kcOHMcVjxBrm4+zvxeQMA3CMpbj1RMLmjN9JxBfsJ2YyU6lGaZVpIm0HncLhHrKz5uNwTE42z+5E97L64qUEd+7CkRMiJm3sM+px62CYwDQaAWAo2ocwu6ioOP3eXYqSmdI2YefMfgDDiL7h9RO7dkNOLsvmzmTs2afIdQwQMQpoA4oMga3QTrDNy7A+SMyIYDVVpLaFnT377HUUZYrStksMoGnWN7wm9++nra6eBdlOBvZ34zZ3YMg82jBw6xruchexniD9Rg8SSUlhDS6XKwXRj3v44URRlCRI64Q9XWR0jPyeLoyG+ehDQ4RiUaxaCGQJ3YaOUwpqy7PRQtAf7yZkl1RVqVtuKNNHRiVsy45daFJSfvFigrt3Y3abAPDKBuzjUxLdtsTjUKSPsF2kfsDp/vsTRVGSIKMStn1nYlXEiy+/lNDOXTjyokgpaDYuwT1+SscR0ZHSIOz3ErdZU5+wR44kiqIkQUYlbKylBW9JGQXuQsZ27SPH4SFCPofJwa1rSItADvgZ0j1oMR1hdVFaWprqsBUlaaaUsEKIm4UQh4UQx4QQ33iTOh8RQhwQQuwXQvxvcsOEgUiM8mNHiDY0YAQCeLsDFJo70I1sWjEo0gWOQhvRLh89shuAwsJKtTi4Mq2cNWGFECbgJ8AtQAOwUgjRcFqdWcA3gSullPOBpE+e3XSkjaKRIUoWLya0bx/B7GJyzf0g82mTOoWGRmWpC8KSbtkFQH3domSHoSgpNZUW9lLgmJSyTUoZBR4H7jitzl8AP5FSDgNIKQeSGyYcGz9+rbskMWHCVJy43C4qZ+CRBnYpqM5JvNZLH4aQ1M2dl+ww3r7FixNFUZJgKv3FCuDEKc+7gMtOqzMbQAixBTAB35FS/un0DQkh7gfuB6iurp5ykHFDEm1uxtA0HA0NeP/zx9gLE3eZa9MvIXd8hDjbkBgYjMYGkHYz1dVpcErn3/891REo08hUWtgz3cNCnvbcDMwClgMrgV8IIfLe8CEpH5ZSNkopG4uKiqYc5B5fkNq2o0Rn1iNMJnzNB8m2DRDHQQuzJuYQm31RRuJebCFJ3GElPz9/yvtQlEwwlYTtAqpOeV4J9JyhzrNSypiU8jhwmEQCJ8UGzwhzO9rIX3wR4YMH8ZndFJo7iInElES3oYFNQ+/z0x3vICdgxp6Xnx73y7nnnkRRlCSYSsLuAGYJIeqEEFbgLuC50+qsAq4FEEK4SXSR25IVZNPBw2SFguRdtIjgzl34sispNHcgDBdtGBTpGkWFdmTEoMPUjdnQKCubcfYNXwhdXYmiKElw1oSVUsaBzwFrgIPAk1LK/UKI7wohbh+vtgYYFEIcADYCX5VSDiYjQE80hti/HwD7wkUEd+2CyhIsWgRdL6dV6hTqguq8xIBTpyXR+M+Zo9YdVqafKZ2klFKuBlaf9toDp/wsgS+Pl6TaOORjXnsr0uHAOnMGoV27sF6RuIHViL4IvwQrgiKrhiF0PHgpxEb9fHVKR5l+0n6m04bBMRZ0tuGcP59YRwcRX4gsUy8SjQOycWLAyRaM45PDaLqOrkFRuboGVpl+0jphdSnZPDDMzM52HIvGj1+zqig0txMzF3KUrMR9dACGQvQGOnAEQc+yItJlZf8rrkgURUmCtJ63t3ssSEHHcUzxGI5FC/Fv3EigZC5uyzNA7vg1sIIclwniknbRRU7AgqM4jW6g/P3vpzoCZRpJk2bozDYMjjGvvRUA+4KFBHfuIlZZR7bJixHPoxUdt6FRmpO40P2Iq5fsoJnyanUrDmV6SuuEXT80xpU9HZgKCsCkEevunlh0LRKdSZvUKdAFpQ4zhsmgxzaIJgWzZl+U4shP8aEPJYqiJEHadok90RhNvhBz24/hWLiQ0O7dxE12nEZilmSPcQnCEFgQZMcN/KYxpEkHzBRXpsGUxJMGk3J2S1GANG5hNw75cIRDZJ3oxL5oIaFduwi463Gb29FNWRykHrcu0ACTP4o31I2IxQDILytPbfCKcp6kbcJuGBzjku5OhJSJEeJduwnNupRCSzuGpShxWw5DI9ckwICOwHHsIYm0mnDk5KY6fEU5L9I2Yf+hvoJvhBLdSWtNDZEjRwgU1lFoPoE08sanJArynYl1nVodHnKCZuzugvSYQ6wo50HaHsOW2CzEjh4mXFVF9PhxkBKDICYRIxwppZU4y3QrxU4zhoxzwuGjetBC8Yw0On4FWLEi1REo00jaJiwk7rLuXLKE4M5dxG1Z2CPHwQHB6Dw6pcFthiBfQMgaZMwWxRUyUV6VZqd0vv3tVEegTCNp2yWOezzEe3uxL1pIcNcuIguupNB8HCnMHJONuAwNGwJLWGcw2kecEAJBYXnV2TeuKBkqbRM21NwCgG3uXMLNzYRqluC2tGM4SjmGA7eukWsWCKDT14mIhgDIL0uzOcS33JIoipIEadslDjU3gckEuo6MxRhzVVEU7UBqVbSh49YFeabE4FJnrA9b2ADS8JROKJTqCJRpJG1bWADnJZcQHm9pg8EwTjFELOqemJJYYDeBS6M3N0JOwIw5y4HNmcL76CjKeZa2LWzxF74AQOf99yNmzcPqOwwFEA1X0Uacm3QzBWZBxBlhxBbHPWghp6QsxVEryvmV1i2s1HVCu/cQnf8e3OZ2AIaji+mTUCwFdl3ijQ3gswXICVgoraxNabyKcr6lbQsLEDl6FMPvJ1A6l8LelzDsBRwN11GgBykwJf6v6fSfIOAcwxmxUlCWhndav+22VEegTCNpnbDBnbsAGNHczLF1gL2cVhKr/E8MOIX6IB4ArOSn4yoTf/M3qY5AmUbSuksc3LUTc2kpw54QedoJ4rp74qL1PLNA5FsYtUrMkTgA+aVpNkKsKEmWtgkrpSS0azfmpZehDR9FI040XEobUdy6Rr5FI5YVZyxLJydgASCvNA0HnZYvTxRFSYK0TdhYVxfxgQFCMxsnBpzCoXpaMaiWAjswEPcwZvWREzTjyM/DYrOnNGZFOd/SNmFPHr/6c+twW9qRJiu98YX4pUa1lrhCp9PfxZh1jNyARU1JVN4V0jZhQ7t3oeXkMBS0Uuw4gcyupRUbheMznKSAzmAffssYeUEbheVpOEKsKEmWtqPEuR/8IM5LLmH7Dh+FluMY2qLxKYka+SaBVmjHNyyJigCWqJU8NeCkvAukbcI6lywhPHsB8T89i614hFCsiFbiuA1BnkWg50niARMiEgHScNL/SR/5SKojUKaRtE1YAM8J38SAUyRYRRtRlhlOrJqgV/fgs/jIDiS+QtpN+j/pr/4q1REo00jaHsMCeDp8FFraAQgE53EcQb1IhNzhP8GYeZScgAWhaeQWl6Qw0rcQDCaKoiRBWifsQIePMtcJZHY5J2QZ0tAo0zSkgHZ/P0HLKLlBCzlFxZjMllSHe2a33pooipIEaZ2wns4x3NZ2DEfNxCV1eWaBLLATMHSCjiCFIScF6Xr8qihJlrYJGw7ECHjHyIqfIG6UjE9JhDyTQC/QkEJjzBbA5Rfkpevxq6Ik2ZQSVghxsxDisBDimBDiG29R78+EEFII0XiugXk6fRSYTyDQx6ckhphrmDALgUd40YVOKOpHi8v0HSFWlCQ7a8IKIUzAT4BbgAZgpRCi4Qz1soHPA9uSEdhAxxhuy3EAIv5a2tCZKxIjwu3+E/jFMNmBxIwnNelfebeYymmdS4FjUso2ACHE48AdwIHT6n0P+FcgKdeTDXYHqMruQlqcjPrqOUGEGqGhC+gY6yOgDU1M+k/rFva++1IdgTKNTKVLXAGcOOV51/hrE4QQS4AqKeXzb7UhIcT9QoidQoidHo/nLXd6wycamF3tQebNoR0LNkNQZBLEc6wE4jHCrjB5QSsmi4Vst3sKXyNF7rtPJa2SNFNJ2DPd90JOvCmEBjwIfOVsG5JSPiylbJRSNhYVFb31TgWYvPvRzZW0oVOia+SYBNE8AQhC2RFKItnklZShjV8MkJa83kRRlCSYSpe4Czj1UphKoOeU59nAAuCl8XvalALPCSFul1LufMeRjXZBeJRYtJhW4lwkTZiEYNg6BsCILUhuMI/8+jQ/fv2zP0s8vvRSSsNQpoeptLA7gFlCiDohhBW4C3ju5JtSylEppVtKWSulrAVeA84tWQH69wMQCZRznCDzSLSixwPtyKgfjz6EZSyuJv0r7ypnTVgpZRz4HLAGOAg8KaXcL4T4rhDi9vMWWX8zAGH/LFoRzEAjCnSM9hEwvDhCGujqlI7y7jKlyf9SytXA6tNee+BN6i4/97AAz2Fkbi1D/Xl4pZ8yk0bYaSKqxwk7AuSk+6R/RTkP0namEx94mMjyx2nFoEBCriYIjy/qH8uX5AatQJqf0lGUJEvfy+s0jVhvmDYMLjbMaGbBiD0AowaB3DjlnnwsdgeuvPxUR/rWPvOZVEegTCPpm7BA9MQQxwmwSDoBaI92YgoH8Fr8zA3lk19anv53W7/zzlRHoEwj6dslBmKeOG1EqRcmQlLS5etGhn30xb04fGTGpP8TJxJFUZIgbVtYGTeI+m20EqZK0/CbBYZhEBHDSF2H0TAFmZCw996beFTnYZUkSNsWNtYfpA+BJs3kaRpBR2JyVbxAkhU0g1SndJR3n/RN2B4fbegsMRKdgGFrAJOhEyk2kzc+QqwmTSjvNumbsO29tBFl8fgMp85YDyIwxli2TlU8Mdk/LW9+pSjnUdomrNk2SIfoZR4aAV3iifSghQL0m0YpjeRgz8rGkZWd6jAV5YJK20GnrIIW2nBSI0yMSInUdETYT29sgJxACfllxakOcWq+ctaLmBRlytI2YSP9RxmSt5CrafRqidtJhm1BdKljGo2SvyhDjl/f975UR6BMI2nbJW5d8HnmyMTg0ogliA2JXuXCpAvio4HMGSE+fDhRFCUJ0raFPRzM4iJpQiLpivWhxcYIVljIPXkv2Ew4Bwvw6U8nHtV5WCUJ0raFPdTnYz4mfAb4GcEYHWYkK0atnjh2zZgWVlGSKG0T1gTMFCZGdEncHMAU9jOgjbx+SidTWlhFSaK0TdhrK/LJFoIRdMxmIBqiO9JHfsiOK78Aq92R6hAV5YJL24T1Hh0CYMgUwmXWiBXZiMs4dp+hWlflXSttB51iJ/wYUjIgBzGHAsRrcgAwBv3kz7woxdG9DX/3d6mOQJlG0jZhKwqc9HWGiFj8GJ4+gjPNWGOCqD+DTukAXH99qiNQppG07RLnNBSxI6gTNweQY8MMOyPMlImucMac0gHYuzdRFCUJ0raFHeoJABJ7loGmx+kTQ8zTS4DBzLq95Be/mHhU52GVJEjbFjbki2JYImTbzRhC0hXupTiSDUKQW1ya6vAUJSXSNmEvvr2CwYIdmCNhqMwlZsTI8pvIcRdjtlpTHZ6ipETaJmxXVxcISWywn1hV4jI6bSSiTuko72ppnbCaphHs7sBfqIGEsGdQJazyrpa2g04Oh4O66ioG9m9nyBGmNlRONBTKrFM6AP/8z6mOQJlG0jZhr7rqKgqMKC8AfQwyiwpgIPPutr5sWaojUKaRtO0SA3hPdCDMJk4Eu6mIFQAZeJXOq68miqIkQdq2sADeznZMdW6iRht5QStjJhM5RRmyNMxJ3/pW4lGdh1WSIK1bWM+JDqKViTtgWccMckvK0ExpfLd1RTnPppSwQoibhRCHhRDHhBDfOMP7XxZCHBBCNAkh1gshas41sLDfj3/Qiz8/ce+c+OAY+aVl57pZRcloZ01YIYQJ+AlwC9AArBRCNJxWbQ/QKKVcBDwF/Ou5BuY90Q7AoD1ImbOMsf7+zDt+VZQkm0oLeylwTErZJqWMAo8Dd5xaQUq5UUoZHH/6GlB5roF5OzsA6DG8zLZUE4+qSROKMpVBpwrg1NuvdQGXvUX9TwEvnOkNIcT9wP0A1dXVb7lT74l2LE4nnYE2ltkbgL7MbGH//d9THYEyjUwlYc90A1Z5xopC3AM0Atec6X0p5cPAwwCNjY1n3MZJM5ZeSqzEQWToQYoiLrxk6L10Fi9OdQTKNDKVLnEXUHXK80qg5/RKQojrgb8FbpdSRs41sBlLL8F+UR0ATp/AbLWRXVB4rpu98NatSxRFSYKptLA7gFlCiDqgG7gL+OipFYQQS4CfAzdLKQeSFVzraGvih+EQ+aVlCC2tz0Kd2T/+Y+JRrTyhJMFZM0BKGQc+B6wBDgJPSin3CyG+K4S4fbzavwFZwO+FEHuFEM8lI7jWkVZKnCX4+vsza5UJRTlPpjTTSUq5Glh92msPnPLzeWk+Wkdaqc+ZyehAH7MuU3NyFSVt+5iGNDg+epwZohxD1zNv0r+inAdpm7Dd/m7CepiyaB6QgZP+FeU8SNvJ/60jiQGn3KCFATL41hw//3mqI1CmkbRPWMtoDJvThSMnN8URvUNz5qQ6AmUaSdsucetIK8XOYgIDXvJKyxHiTPM3MsAf/pAoipIEadvCfn7p5+kP9rP95f+gfPa8VIfzzv3wh4lHdSd2JQnStoUtdZUyP3ceY15P5h6/KkqSpW3CAoz094KUaoRYUcaldcIO9yWmLKuEVZSE9E7Ynm4gg0/pKEqSpe2gE8BIXw/O3DxsTleqQ3nnHn001REo00haJ+xwb09mXgN7qqqqs9dRlClK7y5xX0/md4efeCJRFCUJ0raFjYaCBIaHMn/S/0MPJR7vvDO1cSjTQtq2sMN9vQDkl6sRYkU5KX0Ttnd8hDjTW1hFSaK0TdiR3sQ52Dy1eLiiTEjbhB3u7Sar0I3FZk91KIqSNtJ20Onqez6Jf3go1WGcu6eeSnUEyjSStgnrysvHlZef6jDOndud6giUaSRtu8TTxiOPJIqiJIFK2PNNJaySRCphFSWDqIRVlAyiElZRMohKWEXJIGl7WmfaWL367HUUZYpUwp5vTmeqI1CmEdUlPt9++tNEUZQkUAl7vj35ZKIoShKohFWUDDKlhBVC3CyEOCyEOCaE+MYZ3rcJIZ4Yf3+bEKI22YEqijKFhBVCmICfALcADcBKIUTDadU+BQxLKeuBB4F/SXagiqJMrYW9FDgmpWyTUkaBx4E7TqtzB/Dr8Z+fAlaIjL17laKkr6mc1qkATpzyvAu47M3qSCnjQohRoBDwnlpJCHE/cD9AdXX1Oww5w7z0UqojUKaRqbSwZ2op5Tuog5TyYSllo5SysaioaCrxKYpyiqkkbBdw6mrYlUDPm9URQpiBXGAaLBehKOllKgm7A5glhKgTQliBu4DnTqvzHPDx8Z//DNggpXxDC6soyrk56zHs+DHp54A1gAn4lZRyvxDiu8BOKeVzwC+BR4UQx0i0rHedz6AV5d1qSnOJpZSrgdWnvfbAKT+HgQ8nNzRFUU6nZjopSgZRCasoGUQlrKJkEJWwipJBRKrOvgghPEDHWaq5OW22VIqlWzyQfjGpeN7aVOKpkVKecWZRyhJ2KoQQO6WUjamO46R0iwfSLyYVz1s713hUl1hRMohKWEXJIOmesA+nOoDTpFs8kH4xqXje2jnFk9bHsIqiTJbuLayiKKdQCasoGSRtE/ZsC79d4FiqhBAbhRAHhRD7hRBfSGU8JwkhTEKIPUKI59MgljwhxFNCiEPjv6cr0iCmL43/e7UIIX4nhLBf4P3/SggxIIRoOeW1AiHEWiHE0fHHt3XX8rRM2Cku/HYhxYGvSCnnAZcDn01xPCd9ATiY6iDG/QfwJynlXOAiUhyXEKIC+DzQKKVcQOLS0At92ecjwM2nvfYNYL2Uchawfvz5lKVlwjK1hd8uGCllr5Ry9/jPPhJ/jBWpigdACFEJvBf4RSrjGI8lB7iaxHXRSCmjUsqR1EYFJC4fdYyvguLkjSulnFdSyld448orpy5Y+Gvg/W9nm+masGda+C2lCXLS+JrLS4BtqY2Efwe+BhgpjgNgBuAB/me8i/4LIYQrlQFJKbuBHwCdQC8wKqV8MZUxjSuRUvZCoiEAit/Oh9M1Yae0qNuFJoTIAp4GviilHEthHLcBA1LKXamK4TRmYCnwkJRyCRDgbXb1km382PAOoA4oB1xCiHtSGVMypGvCTmXhtwtKCGEhkay/lVI+k8pYgCuB24UQ7SQOF64TQjyWwni6gC4p5clex1MkEjiVrgeOSyk9UsoY8AywLMUxAfQLIcoAxh8H3s6H0zVhp7Lw2wUzvij6L4GDUsofpSqOk6SU35RSVkop/397d4iCQBBGcfz/qtkD2LyCYBE8h4iYvYDF6kVEEIsXsFtUEIwGNXiLzzATDBq2uDvwfnFg4WOGN7sMw7cd0tzsI6K2t0dEvICHpG4eGgLXuurJ7kBPUiuv35BmHNB9NiwcA7sqDzfy/7C/Gr/VWFIfGAEXSec8Ns+9riyZAau8wd6ASZ3FRMRB0hY4kk75T/z5mqKkNTAA2pKewAJYAhtJU9KmUqkXmq8mmhWkqZ/EZvaFA2tWEAfWrCAOrFlBHFizgjiwZgVxYM0K8gZ78MpBquuXTQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.figure()\n",
    "ax1 = fig.add_subplot(111)\n",
    "sigs = [0, 1, 2, 3, 4, 5, 6,7,8,9,10]\n",
    "for k in num_peaks_all:  \n",
    "    num_peaks = num_peaks_all[k].copy()\n",
    "    num_neurons = np.shape(num_peaks)[0]\n",
    "    ax1.plot(sigs, np.sum(num_peaks[:,:]==1,0)/num_neurons)\n",
    "ax1.set_aspect(1.0/ax1.get_data_ratio())\n",
    "plt.plot([2.75, 2.75], [0,1], c = 'r', ls = '--')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x168ae984888>]"
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAOsAAADrCAYAAACICmHVAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAYN0lEQVR4nO3deXDc9X3/8ed3V7u6LPmUL+SVjTGYyzZGNgWq3yRNJr8kTdIGQkrBQBMSNEnvu7V/82s785OnmWaamWamvzEhaRKb0gZICNM0aUhog1KChWyMscEHttn1LVm27muPb//4rvaQVrJW2u8lvx4zGn+/+/1K37eRX3yvz2GYpomIeF/A7QJEZHoUVhGfUFhFfEJhFfEJhVXEJxRWEZ8oK2bnJUuWmKtXr7apFBHZt2/fJdM06wptKyqsq1evpr29vTRVicgEhmFEJ9umy2ARn1BYRXxCYRXxiaLuWaPRKM3NzQW3bdu2jaamJgBaW1vZs2fPpD9n165dmeWWlhZisVjB/Zqamti2bVvm2Dt37pz0Z27fvp2GhgYA9uzZQ2tra8H9IpEIO3bsyKxP9veBufd3Mk2TBx58iDu23s3gSJJXWl/he8/+C4mkSSKVIpEyreWktfyJ3/lrBkcSDI4meeXZJ+nuOE8ilSKVMhlrUW4Ci67fwMrGD2GaMHj5Aide2p0+XnYfMMGEyAcepnzBMkzgwr4f03PqrYK1li9YRv37H86sn/je30/636lu0wepXbMBgN5TB+k88JMJ+xgYGAasf+CPCRgGQcPg5E92M9J9MbPNAAzDwADq1m3k+rs/QjBgMNB1gcP//k0MrJ0MAwIYhMoMwsEAH3/0t1m9uoHayhB7X/o+Rw+0EQoGCJcFCAcDBANGUb+nyVw1rIZhPAE8ATBv3ryifrgUzzRNUiYkUiku9AwT7+hjcDTJux19nO0eygRpLGDJlMm7gxX0f/cgg6NJBkaS/OzIxfT2bPDGvufFgb1U/HgAgOHThxg4fHHSWs7+14nMcs+FPhK9QwX3G7gyxJVzvQAkevrpGYxP+jNPXx6iLDkIQP9AnJGRZMH9RobinO3OHm9gkv2sHzRCX3rf4f6RKfc92TmQWe7pGyExULjWns4BLpzoAiDRc5me7uFJf+a5V05S9ma/Vcrhk4yczv9vGgwYhIIBak8n2feP/01tZYj5lSFef+8y4ZxQh8qmvtA1iul109jYaOpp8OwMx5Oc7BzgeEcfxy/2c+xiHycvDdA3HGdwNMngaJJkSj2hrlXRL31sn2majYW2FXUZfK1oaWkByLtkKVahUB7v6CfaNYCbWawIBagKl1EZClJdHqQyXEZVKEhVOEhVubVcGba2je2Xu60qHKQ8FMhcLo79GTCsS0lIXyYaY5eWRnodyLncHL997PvsYJqQMk2SKeuqJWWamXXTJP352FfOeir9faaJaZok0+uplMloMkXfcILe4Ti9Qwl6huLp5fTXcILeoTg9Q3ESJfqFK6wFTHa/WYgdoSwLGFZAwmVUlafDEsouV4bK0kFLfx4OTtg2Yb908Mbun8QZpmkyFE9OCHRPTqjHlnuG4nxtip+lsBbhRGc/h872zCqUhgGrFlZx47J5rFtWw7ql81i3tIYlNWGqQmVUhoOEr3LvIv5hGIb1P91wGcvnV1x1/689Nvk2hfUqLg+M8v0DZ3l+/xkOne2d9vcVCuWNy2pYWzePynDQxoplrlJYC0imTM71DPHEt9t5+UjHlPccCqU4RWFNM02Tw+d6eW7fGV44cJaReIrFb+c/gg+XBbh37WLWr6hVKMVx13xYO/tG+P6Bszy37wxHLvQBMBJP5e1zR2QBn7qzno/dvpL5VSE3yhS5NsM6kkjy03c6eG7fGX52rHPCe83yVbeyoDLEE+9by/131rO2To1BxH3XTFhN0+TNMz08v+8ML755jp6hiS1XKkIBPnzrcj71+F9w99rFes0hnjLnw3qhZ5jvvXGW5/ad5kROU7NcW1cv4v47r+Ojt6+gpkKXueJNczKsw/Ek/3H4As/vP8vPj3cWfA9av7CS+zbXc//m62hYXJ23LRq1+v+ONaIX8YI5F9aTnf088vW2vEbgY6rCQT56+wru31zPXWsWEZjkMnesJ0xuTxoRt82psJ6+PMjDT+3lfE9+D4l71i7m/s31fPi25VSXz6m/slxD5sy/3Iu9w2z7ejaolaEgX3jfWu7bfB31C6tcrk5k9uZEWLv6R9j21F6iXVY/yXAwwNcebeSX1y1xuTKR0vF9i/GeoTiPfqON4x1W59+ygME/PrxZQZU5x9dhHRhJ8Jl/auNwepQCw4Cv/MYmPnjLMpcrEyk934Z1OJ7kc99qZ3+sO/PZl+7bwMc3rnSxKhH7+PKedTSR4otP7+cXJ7syn/31x2/h01tWleTnb9++vSQ/R6SUfBfWRDLFH/7rAV4+0pH57M8+fBO/de+akh1DjSHEi3x1GZxKmfz582/xg7fOZz77nfffwBffd4OLVYk4wzdhNU2Tv3rxMM/vP5P57DP3ruaPP3RjyY+1Z8+eKccIFnGDL8JqmiZ/+8Mj7H4tO2fPg1tW8X8/dguGDUPitba2TjqgtohbfBHWr778LrteOZlZ/7VNK2n55O22BFXEqzwf1qdaT/L3Lx3LrH/olmV8+YGN6msq1xxPh/Wf98b4fz94J7PetG4JX33oDkJBT5ctYgvP/qv/3htn2PFCdtKirasX8eQjjZSXaXAyuTZ5Mqw/OnSeP3n2YGYWso318/n6bzVqFEG5pnmuUcR/He3gd595IzOI2frlNXzrs1sdHW4lEok4diyR6fJUWF872UXz7n3Ek1ZQr19Sze7H72JBVdjROmYzIZWIXTxzGfxG7AqPf/N1RhLWmL31Cyt5+vN3UVdT7nJlIt7gibC+fa6Xx77RxsCoNQnustpynv7cXayYX+lyZSLe4XpY3+3o55Gv76V3OAHAouowT3/urgkjDjqpubm56CnkRezmalhPXx5k21N76RoYBaC2oozdj2/lhqU1bpYl4kmuhTWVMvn8t9u50GsNcFYVDvLNz27l1pXz3SpJxNNcC+uxjr7MRFDhsgBPPdbI5shCt8oR8TzXwvr6e1cyyx9Yv5R71mqAM5GpuBfWU5czy42rF7lVhohvuBbW9veyYd2qsIpclSstmM5cGeRceuT86nCQm1d46+nvtm3b3C5BZAJXwtqec7+6uWEhZR7r8tbU1OR2CSITuJKStpxL4C26BBaZFlfCmnu/2rjae69rNAaTeJHjl8FXBkY5djE7L80dq7wX1rGRDXU5LF7i+Jm1PZq9X73tuvnqUC4yTc6HNfeVzRrdr4pMl+NhzX241NjgvUtgEa9yNKxDo0kOne3JrOtJsMj0ORrWA6e7M0O2rFs6j4XVzg7XIuJnjoY1/5WNzqoixXD01U1b3sMl796v7tq1y+0SRCZw7MyaSKbYn/PaprFBZ1aRYjgW1iMX+jIDoq2YX0H9Qg2GJlIMx8LaNq7/qpdngGtpaaGlpcXtMkTyOHbP2h7N7b/q3ftVgFgs5nYJIhM4cmY1TZO2U9n71S1quSRSNEfCGu0a5FL/CGANN3qjhhoVKZojYW0b9341oImQRYrmSFjzB0fz9v2qiFc5EtbcbnEaHE1kZmx/GtzRN8ypSwOANZj37fXeH3Ffnc7Fi2wP676cwdE21S+gvMz7nc01uqF4ke2XwW0eH29JxC9sD2vusKN+eb8ajUaJRqNulyGSx9bL4P6RBIfPWZ3NDQPu9MnIEDt37gTU+0a8xdYz6/7oFVJWX3PWL6+ltiJk5+FE5jRbw5o/n40/zqoiXmVrWMe3XBKRmbMtrKOJFAdOd2fWNTiayOzYFtZD53oYjqcAWLWokuXzK+w6lMg1wbawtmvyKZGSsu3VTW7/Vb+1B96+fbvbJYhMYEtYUykzb2QIvz1camhocLsEkQlsuQw+0dlP92AcgEXVYdbWVdtxGJFrii1hHT+fjZcHRytkz549mWkfRbzClrDmtgf240xxmkxZvMieM+sp/96vinhVycN6rnuIs91DAFSGgty6srbUhxC5JhUX1pG+q+7yes796uaGBYSCjk8BKzInFZekrhMw0j/lLrn3q5rPRqR0ijztmXC2fco9Xs+bKU5hFSmV4q9Ro7+YdFPPYJyjF61L5WDAYNOqBTMuzE2RSIRIJOJ2GSJ5im/BFJs8rPtilzHTnc1vW1lLdbmj07+WzI4dO9wuQWSC4s+sZ16HZLzgptz2wHplI1JaxYc1PgjnDxbcpJ42IvaZ2XuVApfCw/EkB8/0ZNa3+HgYl+bmZpqbm90uQyRPycJ68EwPo0mrs/n1ddUsnlc+q8JEJN/Mwzr2JCkt75WNLoFFSq64sAbST3cHu+DS8bxNag8sYq/iwhrO6ZeacymcTJns10xxIrYqMqzzsss5YT1yoZe+kQQAS2vKWbWosiTFiUhWSc6s4+ez8VtncxE/KK6JUagKyiohMQRX3oPe81C7Im9kiC0+mc9mKpryUbyouDOrYUB9Y3Y99iqmaeY3hpgDjfebmpo0obJ4TvGvbiJ3Z5djr3H68hAXe0cAqCkvY/1ydTYXsUPxLe0bcsIa/QVty3I7my8kGPD//erY+Es6u4qXFB/W+i1gBMBMwcVDvHUiltk0V/qvjo1sqLCKlxR/GVxeA8s3pFdMRk5lnwo3zoGHSyJeNbPmhjn3rfV9bwIQDgbY6NPO5iJ+MLOw5ty3bgkcBeD2+vlUhIIlKUpEJpr1mXWTcYIwcfVfFbHZzMI6byksWgtAuRHnduMkW9foflXETjMe1Dde/0uZ5a3Bo9wZ0ZlVxE4zHtHsveoNrONpAN5X8S7zq0IlK8ptu3btcrsEkQmuemY1DOMJwzDaDcNo7+zszHzeOrIus7zBPAKplD0ViggwjbCapvmkaZqNpmk21tXVZT7/6cUqOs35AFQm+6HzHfuqFJGZ3bPGkyn2x3poS92U/TD6aqlqcl1LSwstLS1ulyGSZ0ZhfftcL0PxJO25YZ1i8G+/icVixGKxq+8o4qAZhXVscLS21Prsh9GJg6iJSOnMKqxHzAjxYHr0iL5z0K2zkYhdig6r1dncGsYlSZDRlbmd0V8rWWEikq/osJ7oHKBrYBSABVUhKtf+cnZjbO48ZBLxmqLDmjuES2PDIgIN+SNHiIg9im7BlDc42uqFcN0tEAhBKg6dR2DwMlT5u+mhOp2LFxUd1vHDjhKugpWbrKkgwXqFs/5XS1agGzS6oXhRUZfB8aRJ7PIgABWhALettFow5Q+iNnfet4p4SVFhHRxNZJY3rVpAuCz97ZH8QdT8LhqNEo1G3S5DJE9Rl8EDIwnC6eW8+Wwi2e5ynD8Ao4PW5bFP7dy5E1DvG/GWos6sA6PJzHLeTHFVi6DuZms5lYCz7SUpTkSyigrrcNwKa8CwxgjOk3t21SsckZKbUXPDW1bWMq983BV0wz3Z5TnUA0fEK2YU1oKDo+WeWc+8DsnExH1EZMZKF9YFEaitt5ZH++HiW7OpS0TGKV1YIf/sOgde4Yh4SdEtmNYsqaauprzwxoa74dBz1nLsF3D3F2dTm2u2b9/udgkiExQd1inns4nkPGSKpTuj+3AW9IaGBrdLEJmg6MvgKSdLrlsPFen5bgY64fLJmdYlIuMUH9appskIBMbdt/rzFc6ePXsy0z6KeEVRYa0Ol7F68VWaEeY1jvDnQ6bW1tbMhMoiXlFUWK+vq8a42j3o+PtWESmJGc91M6mVmyCYflp8+ST0XSz5IUSuRaUPa1k51OcOoqazq0gplD6soM7oIjZQWEV8YsZTPk5p1VYwAmCm4MJbMNwLFbW2HMoOkUjE7RJEJrAnrBW1sOxWK6hmCs60wQ0ftOVQdtixY4fbJYhMYM9lMIx7haPO6CKzZWNY1QNHpJTsC2vuyBFn2yExatuhSq25uZnm5ma3yxDJY19Ya5bDwjXWcmLYGvVQRGbMvrCCXuGIlJC9YW2YW4N/i7jJ2TNrKmXr4UTmMnvDuvgGqFpiLQ93w6Wjth5OZC6zN6yGMSc6o4t4gT0tmHI13ANH/s1ajr0GWx63/ZCzpSkfxYvsD6sPR47QZMriRfZeBgMs3wihamu55zR0n7b9kCJzkf1hDZbBqi3ZdR+0E9YYTOJF9ocVfNc4QqMbihcprCI+4UxY6xshkH6W1fE2DF525LAic4kzYQ1Xw4qN2fXTbY4cVmQucSasMO5SWI0jRIrlUli9/0RYxGscDGtO44iz+yE+5NihReYC+1swjaleAktuhEvHIBW3Arv6XscOX4xdu3a5XYLIBM6dWUH3rSKz4F5Y1RldpCjOhjV35IjTbZBKOnr46WppaaGlpcXtMkTyOBvWBQ1Qs8JaHu2Di4ccPfx0xWIxYrGY22WI5HE2rIahVzgiM+RsWCF/PGGNHCEybc6HNa8z+mtgmo6XIOJHzod16S1QPt9a7r8AV045XoKIHzkf1kDQmhJyjF7hiEyL82GF/Fc4Huzf2tTUpHGYxHOca26Yy+Od0TW6oXiRO2fWlZshGLaWu96F/k5XyhDxE3fCGqqA6+7Mrnvs7BqNRolGo26XIZLHnbDCxFc4HrJz50527tzpdhkieVwMa07jiHd/AokR10oR8QP3wrpqa/a+9dJRePYzkIy7Vo6I17kX1soF8P7t2fWjP4AXvuDZnjgibnMvrAD3/gHc+/vZ9beehX/7AzVBFCnA3bAaBnzwb2DL57Kf7f82/OgvFViRcdwNK1iB/cjfwcaHsp/t/f/wn+r8LZLLnRZM4wUC8ImvQnwQ3n7B+uyVv4NQFTT9kePlbN++/eo7iTjMG2EFa7a5+75mDVF6/D+sz376NxCeB3c94WgpDQ0Njh5PZDrcvwzOVRaGT38L1vyv7Gc//FN4QzO6iXgrrAChSnjwGajP6Ub34u/CoecdK0FTPooXeS+sAOXz4OFnYfkGa91MwXefgKM/dOTwmkxZvMibYQWr0cQjL0Ddems9lYDvPAYn/tPdukRc4t2wAlQvtgK7cI21nhyBf3nIcw3/RZzg7bAC1K6Ax16E2uus9fggPP0AnHvD3bpEHOb9sAIsiMCjL0L1Umt9pBd2fxIuvu1uXSIO8kdYAZbcAI++ABULrPWhK7D716HrhLt1iTjEP2EFWHYrPPJdCNdY6/0X4VufgO7STnURiUSIRCIl/Zkis2WYRTSYb2xsNNvb220sZ5qir8Lu+yCRnpB50fXwmR9CzXJ36xKZJcMw9pmm2Vhom7/OrGMa7oEHn852Xr98Er796zDQ5W5dIjbyZ1gBbvgAPPBNMILWeuc7sOeTMNzjalkidvFvWAHW/yrc9yRgWOvn37Re64wOzOrHNjc309zcPPv6RErI32EFuP1T8Il/yK6f3gvP/CbEh92rScQG/g8rwOZH4cN/m10/9TN49rFZn2FFvGRuhBXgl74Av/J/suvHfgRfvhFe+CK893NIpdyrTaQEvNP5vBSa/sQ6m/78K9b6aD8ceNr6WhCBjb8JGx+0XvWI+MzcObOCNZ7TB/4KPvYVWHJT/rbuGPzsS/APd8A3PmINzDbc606dIjMwt8IKVmAbPwu/vRc+/7I1cuJYE8UxsVetDu1fvhGe/zyceFnjFYvn+bMFU7ESI1bH9TefgeMvgVkgmLXXwYbfgE0P0frOBQDN0SqOm6oF07UR1lz9HXDwO3Dgn6HjcOF96rdY97e33QeVC52tT65pCmshpgkXDsKBZ+Ct78BggaaKwXJY/1HY9DBc/35rBEYRGymsV5MYhXdfss62x35E6ymrQUVTQ0445y2DDZ+Gm38Nlt5sjRMlUmIKazEGumje9knoOs6uXxmcfL/5EVi63hojqm69tbzkJoVYZmWqsOq6brzqxdaZc+nN8IXfs862B78DAx35+/XErK/jP87/XCEWmyisU1l2K/zvFmvyrBMvW7PcnT9gjU5R6IkyKMRiG4V1OoJlcOOHrC+wXgV1nbC65XUcgc7012xCXF1nze0TrrKmDJnucqjKmitI/CmVssYUG+65avdOhXUmysph2S3WV67ZhHhW9VRCuNoKcah6iuUqa73QcqHvKyu3GpnI5EzTGnFzuAeGurOhG85Zznw+/s+edCu66T03UlhLqZQhLkZiyPqa4nnYjBiBwsEuq0iH2LD+NAKFlzOfjS0Hrv59JWdaMzqYKauVmplML6c/y6wns/vlrRfYnkpCKp4NXCphQ90TKaxOuFqILx2zLoVGByE+YHVGKLg8aK2PLcfTX3YxUzDaZ32JfcI11gwUFfOBVyfd7aqvbgzDeAIYm3PxNuBQqWoU2y0BLrldhBTlJtM0awptKOo9q2EY7ZO9AxLv0e/Lf6b6nekxoohPKKwiPlFsWJ+0pQqxi35f/jPp76yoe1YRcY8ug0V8QmEV8QmFVcQnFFYRn1BYRXzifwAMK6VUG/hA2QAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "numneurs = np.zeros((len(sigs), 3))\n",
    "for k in num_peaks_all:    \n",
    "    num_peaks = np.array(num_peaks_all[k]).copy()\n",
    "    for i in sigs:\n",
    "        numneurs[i,0] += np.sum(num_peaks[:,i]==0)\n",
    "        numneurs[i,1] += np.sum(num_peaks[:,i]==1)\n",
    "        numneurs[i,2] += np.sum(num_peaks[:,i]>1)\n",
    "\n",
    "plt.figure()\n",
    "ax = plt.axes()\n",
    "ax.plot(numneurs[:,1:]/np.sum(numneurs[0,:]), lw = 3)\n",
    "ax.set_xlim([0,10])\n",
    "ax.set_ylim([0,1.03])\n",
    "\n",
    "ax.set_xticks([0,5,10])\n",
    "ax.set_xticklabels(['','', ''])\n",
    "ax.set_yticks([0,0.5,1.0])\n",
    "ax.set_yticklabels(['','',''])\n",
    "ax.set_aspect(1.0/ax.get_data_ratio())\n",
    "\n",
    "ax.plot([0,10], [1,1], c = 'k', ls = '--', lw = 2, alpha = 0.6)\n",
    "plt.plot([2.75, 2.75], [0,1], c = 'k', ls = '--', lw = 2, alpha = 0.6)\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
